查看原文
其他

StatQuest生物统计学 - Independent Filtering

冰糖 生信菜鸟团 2022-06-07

目录

什么是Independent Filtering为何进行Independent FilteringedgeR和DESeq2是如何进行Independent Filtering的

什么是Independent Filtering

在edgeR和DESeq2中,原始的read count矩阵除了进行library normalization之外,还会进行Independent Filtering。Independent Filtering筛去一些不太可能有生物学差异的基因,而不管它的p值是否显著,以尽量降低假阴性结果。

为何进行Independent Filtering

对RNA-seq数据进行差异表达分析的时候,如果直接进行假设检验,那么由于基因数目众多,因而进行的多次统计学检验会增加假阳性的风险。以一般的统计学显著cutoff值0.05来说,每一次的假设检验都会有 5%的假阳性,随着假设检验的增加,那么假阳性风险会非常高,这种情形叫做多重比较。

当然实际情形下并不会出现很高的假阳性概率,原因在于edgeR和DESeq2都会对假设检验的p值进行FDR校正。所以假阳性的风险可以维持在5%以下。

FDR是一种校正思想,而不是具体的方法,下一周将会介绍FDR思想和其的一种具体实现方法Benjamini-Hochberg方法。

然而,虽然FDR可以控制假阳性概率,但是随着基因数目的增加,大量真实差异基因会被掩盖,也就是说假阴性概率提高了,而FDR是无法修正这种情形的。

而Independent Filtering就是为了降低假设检验的假阴性。

下面的例子说明了假阴性为何会随着基因数目的增加而增加。

假如从两个不同的正态总体中随机各挑选3个数字,并进行假设检验,差异显著性为0.05,那么会有5%的假阳性。进行上述模拟抽样统计,结果如下,949次抽样是显著的,而有51次抽样是得出两个样本不显著(I类统计学错误)。

上述的1000个抽样由于来自不同的正态总体,所以他们实际是有差异的,代表着差异基因。然后继续进行1000次模拟抽样,这次抽样来自同一总体,由于是同一总体,因此是代表无差异基因。根据5%的假阳性,那么会有约50次抽样显著,即便他们是来自于同一总体,原本是没有差异的。

下图代表合并上述两次抽样结果,1000个来自于不同总体的抽样和1000次来自于同一总体的抽样

可以看出此时的阳性有抽样共有993个,其中真阳性是949个(来自于不同总体),假阳性是44个(来自于同一总体)。

由于是多重比较,实施FDR之后,则有846个阳性结果保留下来,其中真阳性827个,假阳性19个,假阳性率2%。

继续进行进行试验,将1000个来自于同一总体的抽样次数增加为5000个,下图为为其结果,1000个来自于不同总体的抽样和5000次来自于同一总体的抽样,并且已实施FDR:

此时的阳性结果256个,其中真阳性250个,假阳性6个,假阳性率2%。

继续进行进行试验,将5000个来自于同一总体的抽样次数增加为10,000个,下图为为其结果,1000个来自于不同总体的抽样和10,000次来自于同一总体的抽样,并且已实施FDR:

此时的阳性结果56个,其中真阳性54个,假阳性2个,假阳性率2%。

可以发现,随着模拟无差异基因的假设检验的数目的增加(1000, 5000, 10000),假阳性可以控制的非常好:2%,2%,4%,但是真阳性结果却被大量掩盖,从89%,26%到6%。所以从这个模拟实验中,可以非常合理的看出在RNA-seq分析中,如果不进行任何处理的话,将会有大量的真实差异表达基因被掩盖,假阴性增加。

edgeR和DESeq2是如何进行Independent Filtering的

幸运的是,Independent Filtering可以降低这种假阴性结果。

Independent Filtering的基本思想是将超低Read数的基因从数据集中移除,因为过低的Read数(如只有1、2个Read数)首先不太可能有生物学意义,其次即便有生物学意义,也不太可能测得准。

edgeR的Independent Filtering方法

edgeR会将在2个及以上样本中的CPM值大于1的基因保留,而剔除剩余的基因。

对于以下已转为CPM的read count矩阵,gene1在两个样本中CPM大于1,因此保留,gene2同理。gene3虽然在一个样本中有一个极大的CPM值103.3,但是由于只有一个样本CPM大于1,因此剔除,gene4、5同理。gene6在两个样本中CPM值大于1,因此也会保留,即使剩余的两个样本中CPM值为0。

然而需要说明的是,将CPM的cutoff值定位1,是非常武断的做法,这个CPM值会因到测序深度而有很大的影响,因此尝试几个CPM值是更好的方案。

什么是CPM?

CPM: Counts Per Million,它其实是只进行RPKM的第一部分,只进行测序深度标准化。

对于每一个Read,除以这个Read所在样本的总Read数,再乘以10^6即可。

举例来说,有以下Read count数据:

Control

gene1 65

gene2 5

gene3 0

... ...

样本Control共有5324112个Read。

那么gene1的CPM就是65/5324112*10^6=12.3,gene2的CPM是5/5324112*10^6=0.9,其他类推。

DESeq2的Independent Filtering方法

DESeq2不同于edgeR,它是根据每一个基因的平均CPM值决定的。

而且它会自动尝试众多的cutoff值,并在这些cutoff值给出最佳的cutoff值。具体来说,就是DESeq2会将根据CPM值将gene排序,然后尝试不同的分位数,将次分位数一下的基因剔除,然后计算此时的显著基因数量,并对这个“显著基因数量-分位数”曲线进行拟合(见下图,y为显著基因数量,下x轴为分位数,上x轴为分位数对应的CPM cutoff值),选出最佳的分位数。

参考资料

StatQuest课程:https://statquest.org/video-index/


猜你喜欢

生信基础知识100讲

生信菜鸟团-专题学习目录(5)

生信菜鸟团-专题学习目录(6)

生信菜鸟团-专题学习目录(7)

还有更多文章,请移步公众号阅读

▼ 如果你生信基本技能已经入门,需要提高自己,请关注下面的生信技能树,看我们是如何完善生信技能,成为一个生信全栈工程师。

▼ 如果你是初学者,请关注下面的生信菜鸟团,了解生信基础名词,概念,扎实的打好基础,争取早日入门。



      

    


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存